
global gammagrid elementsize lowerbound table m n agent Y rho beta P S
agent=4;
lowerbound=0.05;
elementsize=0.05;
upperbound=(1-(agent-1)*lowerbound);

theta=[lowerbound:elementsize:upperbound];
m=length(theta);
for r1=1:m
index4=0;
for k=0:r1-2
index4=index4+sum(1:m-k);
end
for r2=1:m-(r1-1)
index3=sum(m-(r1-1):-1:(m-(r1-1))-(r2-2)); 
for r3=1:m-(r1-1)-(r2-1)
table(r1,r2,r3)=index3+index4+r3;
end
end
end
size(table)
for r1=1:m
for r2=1:m-(r1-2)
for r3=1:m-(r1-1)-(r2-1)
gammagrid(table(r1,r2,r3),1:3)=[(r1-1)*elementsize+lowerbound (r2-1)*elementsize+lowerbound (r3-1)*elementsize+lowerbound ];
end
end
end
gammagrid(:,4)=1-sum(gammagrid(:,1:3),2);
gammagrid(:,5)=1;
gammagrid(:,6)=gammagrid(:,2)./gammagrid(:,1);
gammagrid(:,7)=gammagrid(:,3)./gammagrid(:,1);
gammagrid(:,8)=gammagrid(:,4)./gammagrid(:,1);

n=length(gammagrid);

d1=n;   %length of gammagrid
d2=16;  %states
d3=agent;   %agent





S=gridmake(incomevals,incomevals,incomevals,incomevals);
SS=gridmake([1:length(incomevals(:,1))]', [1:length(incomevals(:,1))]', [1:length(incomevals(:,1))]',[1:length(incomevals(:,1))]');
state=length(S);
Y=sum(S,2);				% possible aggregate income outcome


for k=1:state
    for j=1:state
        P(j,k)=Q1(SS(j,1),SS(k,1))*Q1(SS(j,2),SS(k,2))*Q1(SS(j,3),SS(k,3))*Q1(SS(j,4),SS(k,4));
        %P(j,k)=0.0625;
    end
end

% possible aggregate income outcome



w=zeros(n,state,agent);
c=zeros(n,state,agent);
caut=zeros(state,agent);
%%%Compute the set of autarky values
for k=1:agent
    caut(:,k)=S(:,k);
end
waut=zeros(state,agent);
wautnew=zeros(state,agent);
dwaut=1;
crit=0.0000000001;
while dwaut>crit
    
    wautnew(:,1)=utility(rho,caut(:,1))+beta*P*waut(:,1);
    wautnew(:,2)=utility(rho,caut(:,2))+beta*P*waut(:,2);
    wautnew(:,3)=utility(rho,caut(:,3))+beta*P*waut(:,3);
    wautnew(:,4)=utility(rho,caut(:,4))+beta*P*waut(:,4);
    dwaut=max(max(abs(wautnew-waut)));
    waut=wautnew;
end

% Note that phi contains the Lagrange multipliers on the enforcement
% constraint
phi=zeros(n,state,agent);

for j=1:state
    for g=1:agent
c(:,j,g)=Y(j)*gammagrid(:,g);
     end
end

cfirstbest=c;

% Compute the value functions
w=zeros(n,state,agent);
wnew=w;
dw=1;
while dw>crit
    for j=1:state
        wnew(:,j,1)=utility(rho,cfirstbest(:,j,1))+beta*w(:,:,1)*P(j,:)';
        wnew(:,j,2)=utility(rho,cfirstbest(:,j,2))+beta*w(:,:,2)*P(j,:)';
        wnew(:,j,3)=utility(rho,cfirstbest(:,j,3))+beta*w(:,:,3)*P(j,:)';
        wnew(:,j,4)=utility(rho,cfirstbest(:,j,4))+beta*w(:,:,4)*P(j,:)';

    end
            dw=max(max(max(abs(w-wnew))));
        w=wnew;
end
    wfirstbest=w;
    ewfirstbest=zeros(n,state,agent);
    %%%
    for j=1:state
        for g=1:agent
            ewfirstbest(:,j,g)=wfirstbest(:,:,g)*P(j,:)';
        end
    end
    ewsb=ewfirstbest;
    wsb=wfirstbest;


gap=1;
gapcritsetexact=0.1
gapend=0.001;
options=optimset('Display','iter');
scalefactor=1;

count=0;
eps=0.0000001;
t=0;
xkeep=[];
outputkeep=[];
xaux=[];
csol=[];
phisol=[];
wsbnewsol=[];
wsbnew=wsb;


for t=1:20


agentperm=perms(1:agent);
count=0
countgrid=ones(n,state);
countperm=ones(state,1);


for j=1:state
stateperm=perms(S(j,:));

index=[1 j]    
%%%%calculate the allocations in each state when $n-1$ agents have a binding constraint    
rr=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)<=waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)<=waut(j,2) & utility(rho,cfirstbest(:,j,3))+beta*ewsb(:,j,3)<=waut(j,3));
if size(rr)>0 
x0=[cfirstbest(rr(end),j,1) cfirstbest(rr(end),j,2) cfirstbest(rr(end),j,3)];
            bind=[1 2 3];
            nobind=4;
            [xsol,output]=agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
                    
            outputkeepsol(j,1:agent)=output(1:4);
            xkeepsol(j,1:agent)=[xsol(1:3) Y(j)-xsol(1)-xsol(2)-xsol(3)];
                        
            if output(6)<0
                error('did not converge')
            end
for i=1:length(rr)
index=[rr(i),j];
if countgrid(rr(i),j)>0

   
gridperm=perms(gammagrid(rr(i),1:4));
for k=1:length(gridperm)
iperm(k)=find(gammagrid(:,1)<=gridperm(k,1)+eps & gammagrid(:,1)>=gridperm(k,1)-eps& gammagrid(:,2)<=gridperm(k,2)+eps & gammagrid(:,2)>=gridperm(k,2)-eps & gammagrid(:,3)<=gridperm(k,3)+eps & gammagrid(:,3)>=gridperm(k,3)-eps & gammagrid(:,4)<=gridperm(k,4)+eps & gammagrid(:,4)>=gridperm(k,4)-eps);
jperm(k)=find(S(:,1)==stateperm(k,1) & S(:,2)==stateperm(k,2) & S(:,3)==stateperm(k,3) & S(:,4)==stateperm(k,4));
end
            
csol(rr(i),j,1:3)=xkeepsol(j,1:3,1);
csol(rr(i),j,4)=Y(j)-xkeepsol(j,1,1)-xkeepsol(j,2,1)-xkeepsol(j,3,1);
wsbnewsol(rr(i),j,1:agent)=outputkeepsol(j,1:agent,1);
countgrid(rr(i),j)=0;

for g=1:agent
gammar(rr(i),j,g)=csol(rr(i),j,g)/csol(rr(i),j,1);
end
for g=1:agent
phisol(rr(i),j,g)=((gammagrid(rr(i),nobind(end))/gammagrid(rr(i),g))/(gammar(rr(i),j,nobind(end))/gammar(rr(i),j,g)))^rho-1;
end

%identical cases

for k=1:length(gridperm)
    for g=1:agent
    c(iperm(k),jperm(k),g)=csol(rr(i),j,agentperm(k,g));
    phi(iperm(k),jperm(k),g)=phisol(rr(i),j,agentperm(k,g));
    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(rr(i),j,agentperm(k,g));
    end
    for g=1:agent
    gammar(iperm(k),jperm(k),g)=c(iperm(k),jperm(k),g)/c(iperm(k),jperm(k),1);    
    end
countgrid(iperm(k),jperm(k))=0;
end

end
end


%%%%and update the xkeeps for below
if countperm(j)>0
agent3perm=perms(1:agent-1);
state3perm=perms(S(j,1:agent-1));
for k=1:length(agent3perm)
j3perm(k)=find(S(:,1)==state3perm(k,1) & S(:,2)==state3perm(k,2) & S(:,3)==state3perm(k,3) & S(:,4)==S(j,4));
end
for k=1:length(agent3perm)
for g=1:agent-1
xkeep(j3perm(k),g)=xkeepsol(j,agent3perm(k,g));
outputkeep(j3perm(k),g)=outputkeepsol(j,agent3perm(k,g));
end
xkeep(j3perm,4)=xkeepsol(j,4);
outputkeep(j3perm,4)=outputkeepsol(j,4);
countperm(j3perm(k))=0;
end
end % end of countperm if loop

end

end % end of j loop

clear outputkeepsol xkeepsol
%%%%allocate the other permutations
BB=nchoosek(1:4,3);
for k=1:size(BB,1)
   NB(k,:)=find(ismember([1:4],BB(k,:))==0); 
end
B=[BB NB];
for j=1:state
    for k=2:size(BB,1)
    jjperm=find(S(:,B(k,1))==S(j,B(1,1)) & S(:,B(k,2))==S(j,B(1,2)) & S(:,B(k,3))==S(j,B(1,3))  & S(:,B(k,4))==S(j,B(1,4)));
    for g=1:agent
    outputkeep(jjperm,B(k,g),k)=outputkeep(j,B(1,g),1);
    xkeep(jjperm,B(k,g),k)=xkeep(j,B(1,g),1);
    end
    end
end



%for j=1:state
tic
for j=1:state
stateperm=perms(S(j,:));
count=count+1

for i=1:n
 
index=[i,j];
if countgrid(i,j)>0

   
gridperm=perms(gammagrid(i,1:4));
for k=1:length(gridperm)
iperm(k)=find(gammagrid(:,1)<=gridperm(k,1)+eps & gammagrid(:,1)>=gridperm(k,1)-eps& gammagrid(:,2)<=gridperm(k,2)+eps & gammagrid(:,2)>=gridperm(k,2)-eps & gammagrid(:,3)<=gridperm(k,3)+eps & gammagrid(:,3)>=gridperm(k,3)-eps & gammagrid(:,4)<=gridperm(k,4)+eps & gammagrid(:,4)>=gridperm(k,4)-eps);
jperm(k)=find(S(:,1)==stateperm(k,1) & S(:,2)==stateperm(k,2) & S(:,3)==stateperm(k,3) & S(:,4)==stateperm(k,4));
end

% Constraint binding for
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%   Person 1 only  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1))<(utility(rho,S(j,1))+beta*P(j,:)*waut(:,1)) ...
    && (utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2))>=(utility(rho,S(j,2))+beta*P(j,:)*waut(:,2)) ...
    && (utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3))>=(utility(rho,S(j,3))+beta*P(j,:)*waut(:,3))...
    && (utility(rho,cfirstbest(i,j,4))+beta*ewsb(i,j,4))>=(utility(rho,S(j,4))+beta*P(j,:)*waut(:,4));
bind=1;
nobind=[2 3 4];
XBIND=[xkeep(j,1:3,1) xkeep(j,4,2)];
clear x0
xaux(1) = XBIND(1);
xaux(2) = (Y(j)-XBIND(1))/(1+gammagrid(i,7)/gammagrid(i,6)+gammagrid(i,8)/gammagrid(i,6));
xaux(3) = gammagrid(i,7)/gammagrid(i,6)*xaux(2);
xaux(4) = gammagrid(i,8)/gammagrid(i,6)*xaux(2);

%%%%%now go through all the possibilities
bindpost=bind;
nobindpost=[];
for g=1:length(nobind)
     if xaux(nobind(g))<XBIND(nobind(g))
     bindpost=[bindpost nobind(g)];
     else
     nobindpost=[nobindpost nobind(g)];
     end
end


%%%%%now go through the possibilities    
if length(bind)==length(bindpost)

index=[i j];
x0=[xaux(1:3)];
[x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
if output(6)<=0
%%%this is on case by case basis
x0=xkeep(8,1:3,2);
bindaux=[1 4]
nobindaux=[2 3]
[x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bindaux,nobindaux);
x0=x;
[x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
end
if output(6)<=0
    if i==106 & j==8 & (beta==0.86 | beta==0.94 | beta==0.84)
    rr=find(gammagrid(:,1)>=0.05 & gammagrid(:,2)>=0.35 & gammagrid(:,3)>=0.35);
    x0=[squeeze(c(rr(1),j,1:agent-1))']
    [x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
    end
    if i==371 & j==8 &( beta==0.86 | beta==0.96)
    rr=find(gammagrid(:,1)>=0.10 & gammagrid(:,2)>=0.35 & gammagrid(:,3)>=0.35);
    x0=[squeeze(c(rr(1),j,1:agent-1))']
    [x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
    end
    if output(6)<=0
    x0=squeeze(c(i-1,j,1:3));
    [x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
    end
end
if output(6)<=0
error('Did not converge');
end

elseif length(bind)<length(bindpost)
    %%%%now checking for any two who are constrained 
       
    if length(bindpost)==2
        for g=1:length(bindpost)
        xaux(bindpost(g)) = XBIND(bindpost(g));
        end
        for g=1:length(nobindpost)
        xaux(nobindpost(1)) = (Y(j)-XBIND(bindpost(1))-XBIND(bindpost(2)))/(1+gammagrid(i,nobindpost(2))/gammagrid(i,nobindpost(1)));
        xaux(nobindpost(2)) = gammagrid(i,nobindpost(2))/gammagrid(i,nobindpost(1))*xaux(nobindpost(1));
        end
        %%%%now check if none of the other non-binding ones has a constraint.
        nobindpost2=[];
        bindpost2=bindpost;
        
        for g=1:length(nobindpost)
            if xaux(nobindpost(g))<XBIND(nobindpost(g))
            bindpost2=[bindpost2 nobindpost(g)];
            
            else
            nobindpost2=[nobindpost2 nobindpost(g)];
            end
        end
        bindpost2=sort(bindpost2);
        nobindpost2=sort(nobindpost2);
      %%%%if it's incentive compatible for everyone, then solve
        if length(bindpost)==length(bindpost2)
        bind=bindpost2;
        nobind=nobindpost2;
      
        
        rr=find(bind==4);
        if size(rr)==0
        index=[i j];
        x0=[xaux(1:3)];
        [x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
        
        if output(6)<=0 & model>1
        x0=[CKEEP(i,j,1,t,model-1) CKEEP(i,j,2,t,model-1) CKEEP(i,j,3,t,model-1)] ;
        [x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);

        end
        if output(6)<=0
        error('Did not converge');
        end
        end
   
        elseif length(bindpost)<length(bindpost2)
        bind=bindpost2;
        nobind=nobindpost2;
        rr=find([bindpost2(1)==BB(:,1) & bindpost2(2)==BB(:,2) & bindpost2(3)==BB(:,3)]);  
        x=xkeep(j,:,rr); 
        output=outputkeep(j,:,rr);
            
        end
    elseif length(bindpost)==3
        bind=bindpost;
        nobind=nobindpost;
        rr=find([bindpost(1)==BB(:,1) & bindpost(2)==BB(:,2) & bindpost(3)==BB(:,3)]);  
        x=xkeep(j,:,rr); 
        output=outputkeep(j,:,rr);
    
        
        
    end
    
end

        
 

rr=find(bind==4);
if size(rr)==1 
    
else

csol(i,j,1:agent-1)=x(1:agent-1);
csol(i,j,agent)=Y(j)-sum(x(1:agent-1));
wsbnewsol(i,j,1:agent)=output(1:4);

for g=1:agent
gammar(i,j,g)=csol(i,j,g)/csol(i,j,1);
end
for g=1:agent
phisol(i,j,g)=((gammagrid(i,nobind(end))/gammagrid(i,g))/(gammar(i,j,nobind(end))/gammar(i,j,g)))^rho-1;
end

%identical cases

for k=1:length(gridperm)
    for g=1:agent
    c(iperm(k),jperm(k),g)=csol(i,j,agentperm(k,g));
    phi(iperm(k),jperm(k),g)=phisol(i,j,agentperm(k,g));
    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(i,j,agentperm(k,g));
    end
    for g=1:agent
    gammar(iperm(k),jperm(k),g)=c(iperm(k),jperm(k),g)/c(iperm(k),jperm(k),1);    
    end
countgrid(iperm(k),jperm(k))=0;
end

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Constraint binding for Person 1 and Person 2 %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif  scalefactor*(utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1))<scalefactor*(utility(rho,S(j,1))+beta*P(j,:)*waut(:,1)) ...
    && scalefactor*(utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2))<scalefactor*(utility(rho,S(j,2))+beta*P(j,:)*waut(:,2)) ...
    && scalefactor*(utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3))>=scalefactor*(utility(rho,S(j,3))+beta*P(j,:)*waut(:,3))...
    && scalefactor*(utility(rho,cfirstbest(i,j,4))+beta*ewsb(i,j,4))>=scalefactor*(utility(rho,S(j,4))+beta*P(j,:)*waut(:,4));


bind  =[ 1 2];
nobind=[ 3 4];
XBIND=[xkeep(j,1:3,1) xkeep(j,4,2)];
clear x0
xaux(1) = XBIND(1);
xaux(2) = XBIND(2);
xaux(3) = (Y(j)-XBIND(1)-XBIND(2))/(1+gammagrid(i,8)/gammagrid(i,7));
xaux(4) = gammagrid(i,8)/gammagrid(i,7)*xaux(3);

%%%%%now go through all the possibilities
bindpost=bind;
nobindpost=[];
for g=1:length(nobind)
     if xaux(nobind(g))<XBIND(nobind(g))
     bindpost=[bindpost nobind(g)];
     else
     nobindpost=[nobindpost nobind(g)];
     end
end
  
if length(bind)==length(bindpost)

index=[i j];
x0=[xaux(1:3)];
[x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
if output(6)<=0
   %%%this is on case by case basis
bindaux=[1 2 4];
nobindaux=3;
[x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bindaux,nobindaux);
x0=x;
[x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
if output(6)<=0
    if i==59 & j==8 & beta==0.86
    rr=find(gammagrid(:,1)>=0.05 & gammagrid(:,2)>=0.15 & gammagrid(:,3)>=0.6);
    x0=[squeeze(c(rr(1),j,1:agent-1))']
    [x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
    end
    if i==45 & j==8 & beta==0.84
    rr=find(gammagrid(:,1)>=0.05 & gammagrid(:,2)>=0.1 & gammagrid(:,3)>=0.6);
    x0=[squeeze(c(rr(1),j,1:agent-1))']
    [x, output] = agentmaxrho_4agents_V2(x0,index,ewsb,waut,bind,nobind);
    end
    if i==59 & j==8 & beta==0.84
    rr=find(gammagrid(:,1)>=0.05 & gammagrid(:,2)>=0.15 & gammagrid(:,3)>=0.6);
    x0=[squeeze(c(rr(1),j,1:agent-1))']
    [x, output] = agentmaxrho_4agents_V2(x0,index,ewsb,waut,bind,nobind);
    end
    if i==181 & j==8 & beta==0.84
    rr=find(gammagrid(:,1)>=0.1 & gammagrid(:,2)>=0.1 & gammagrid(:,3)>=0.55);
    x0=[squeeze(c(rr(1),j,1:agent-1))']
    [x, output] = agentmaxrho_4agents_V2(x0,index,ewsb,waut,bind,nobind);
    end
    if i==221 & j==8 & beta==0.84
    rr=find(gammagrid(:,1)>=0.1 & gammagrid(:,2)>=0.2 & gammagrid(:,4)>=0.15);
    x0=[squeeze(c(rr(1),j,1:agent-1))']
    [x, output] = agentmaxrho_4agents_V2(x0,index,ewsb,waut,bind,nobind);
    end
    if i==352 & j==8 & beta==0.84
    rr=find(gammagrid(:,1)>=0.15 & gammagrid(:,2)>=0.2 & gammagrid(:,4)>=0.15);
    x0=[squeeze(c(rr(1),j,1:agent-1))']
    [x, output] = agentmaxrho_4agents_V2(x0,index,ewsb,waut,bind,nobind);
    end
    if i==457 & j==8 & beta==0.84
    rr=find(gammagrid(:,1)>=0.15 & gammagrid(:,2)>=0.2 & gammagrid(:,4)>=0.15);
    x0=[squeeze(c(rr(1),j,1:agent-1))']
    [x, output] = agentmaxrho_4agents_V2(x0,index,ewsb,waut,bind,nobind);
    end
    if i==13 & j==8 & beta==0.96
    x0=xkeep(j,1:3,1)
    [x, output] = agentmaxrho_4agents_V2(x0,index,ewsb,waut,bind,nobind);
    end
    if i==14 & j==8 & beta==0.96
    x0=squeeze(c(i-1,j,1:3));
    [x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
    end
    if i==29 & j==8 & beta==0.96
    x0=squeeze(c(13,j,1:3));
    [x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
    end
    if i==30 & j==8 & beta==0.96
    x0=squeeze(c(i-1,j,1:3));
    [x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
    end
    if i==44 & j==8 & beta==0.96
    x0=squeeze(c(i-1,j,1:3));
    [x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
    end
    if i==195 & j==8 & (beta==0.86 | beta==0.84)
    rr=find(gammagrid(:,1)>=0.1 & gammagrid(:,2)>=0.1 & gammagrid(:,3)>=0.55);
    x0=[squeeze(c(rr(1),j,1:agent-1))']
    [x, output] = agentmaxrho_4agents_V2(x0,index,ewsb,waut,bind,nobind);
    end
    if i==467 & j==8 & ( beta==0.96 | beta==0.84)
    rr=find(gammagrid(:,1)>=0.2 & gammagrid(:,2)>=0.2 & gammagrid(:,3)>=0.4);
    x0=[squeeze(c(rr(1),j,1:agent-1))']
    [x, output] = agentmaxrho_4agents_V2(x0,index,ewsb,waut,bind,nobind);
    end
    if output(6)<=0
    x0=squeeze(c(i-1,j,1:3));
    [x, output] = agentsolverho_4agents(x0,index,ewsb,waut,bind,nobind);
    end

end
end
   
if output(6)<=0
error('Did not converge');
end

elseif length(bind)<length(bindpost)
    %%%%that means three people are constrained
  
        bind=bindpost;
        nobind=nobindpost;
        rr=find([bindpost(1)==BB(:,1) & bindpost(2)==BB(:,2) & bindpost(3)==BB(:,3)]);  
        x=xkeep(j,:,rr); 
        output=outputkeep(j,:,rr);
       
   
end




csol(i,j,1:agent-1)=x(1:agent-1);
csol(i,j,agent)=Y(j)-sum(x(1:agent-1));
wsbnewsol(i,j,1:agent)=output(1:4);

for g=1:agent
gammar(i,j,g)=csol(i,j,g)/csol(i,j,1);
end
for g=1:agent
phisol(i,j,g)=((gammagrid(i,nobind(end))/gammagrid(i,g))/(gammar(i,j,nobind(end))/gammar(i,j,g)))^rho-1;
end

%identical cases

for k=1:length(gridperm)
    for g=1:agent
    c(iperm(k),jperm(k),g)=csol(i,j,agentperm(k,g));
    phi(iperm(k),jperm(k),g)=phisol(i,j,agentperm(k,g));
    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(i,j,agentperm(k,g));
    end
    for g=1:agent
    gammar(iperm(k),jperm(k),g)=c(iperm(k),jperm(k),g)/c(iperm(k),jperm(k),1);    
    end
countgrid(iperm(k),jperm(k))=0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constraint binding for Person 1 and Person 2 and Person 3%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif scalefactor*(utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1))<scalefactor*(utility(rho,S(j,1))+beta*P(j,:)*waut(:,1)) ...
    && scalefactor*(utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2))<scalefactor*(utility(rho,S(j,2))+beta*P(j,:)*waut(:,2)) ...
    && scalefactor*(utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3))<scalefactor*(utility(rho,S(j,3))+beta*P(j,:)*waut(:,3))...
    && scalefactor*(utility(rho,cfirstbest(i,j,4))+beta*ewsb(i,j,4))>=scalefactor*(utility(rho,S(j,4))+beta*P(j,:)*waut(:,4));

bind=[1 2 3];
nobind=4;
csol(i,j,1:3)=xkeep(j,1:3,1);
csol(i,j,4)=Y(j)-xkeep(j,1,1)-xkeep(j,2,1)-xkeep(j,3,1);
wsbnewsol(i,j,1:agent)=outputkeep(j,1:agent,1);
countgrid(i,j)=0;

for g=1:agent
gammar(i,j,g)=csol(i,j,g)/csol(i,j,1);
end
for g=1:agent
phisol(i,j,g)=((gammagrid(i,nobind(end))/gammagrid(i,g))/(gammar(i,j,nobind(end))/gammar(i,j,g)))^rho-1;
end

%identical cases

for k=1:length(gridperm)
    for g=1:agent
    c(iperm(k),jperm(k),g)=csol(i,j,agentperm(k,g));
    phi(iperm(k),jperm(k),g)=phisol(i,j,agentperm(k,g));
    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(i,j,agentperm(k,g));
    end
    for g=1:agent
    gammar(iperm(k),jperm(k),g)=c(iperm(k),jperm(k),g)/c(iperm(k),jperm(k),1);    
    end
countgrid(iperm(k),jperm(k))=0;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%Constraint binding for noone%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif scalefactor*(utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1))>=scalefactor*(utility(rho,S(j,1))+beta*P(j,:)*waut(:,1)) ...
    && scalefactor*(utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2))>=scalefactor*(utility(rho,S(j,2))+beta*P(j,:)*waut(:,2)) ...
    && scalefactor*(utility(rho,cfirstbest(i,j,3))+beta*ewsb(i,j,3))>=scalefactor*(utility(rho,S(j,3))+beta*P(j,:)*waut(:,3))...
    && scalefactor*(utility(rho,cfirstbest(i,j,4))+beta*ewsb(i,j,4))>=scalefactor*(utility(rho,S(j,4))+beta*P(j,:)*waut(:,4));
bind=[];
nobind=[1:4];
csol(i,j,:)=cfirstbest(i,j,:);
for g=1:agent
wsbnewsol(i,j,g)=utility(rho,cfirstbest(i,j,g))+beta*ewsb(i,j,g);
end
countgrid(i,j)=0;

for g=1:agent
gammar(i,j,g)=csol(i,j,g)/csol(i,j,1);
end
for g=1:agent
phisol(i,j,g)=((gammagrid(i,nobind(end))/gammagrid(i,g))/(gammar(i,j,nobind(end))/gammar(i,j,g)))^rho-1;
end

%identical cases

for k=1:length(gridperm)
    for g=1:agent
    c(iperm(k),jperm(k),g)=csol(i,j,agentperm(k,g));
    phi(iperm(k),jperm(k),g)=phisol(i,j,agentperm(k,g));
    wsbnew(iperm(k),jperm(k),g)=wsbnewsol(i,j,agentperm(k,g));
    end
    for g=1:agent
    gammar(iperm(k),jperm(k),g)=c(iperm(k),jperm(k),g)/c(iperm(k),jperm(k),1);    
    end
countgrid(iperm(k),jperm(k))=0;
end


end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%All the updates%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





end % end of countgrid if loop
end  % end of i loop


end  % end of j loop
toc
gap=scalefactor*norm((wsbnew(:,:,1)-wsb(:,:,1))./wsb(:,:,1));
gap
%gamma's
for j=1:state
for g=1:agent
gammar(:,j,g)=c(:,j,g)./c(:,j,1);
end
end
wsb=wsbnew;
for j=1:state
        for g=1:agent
ewsb(:,j,g)=wsb(:,:,g)*P(j,:)';
        end
end
gap
if sum(sum(countgrid))>0
    error('countgrid not zero')
end

end % of iteration loop


